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ABSTRACT 



The first objects to arise in a cold dark matter universe present a daunting 
challenge for models of structure formation. In the ultra small-scale limit, CDM 
structures form nearly simultaneously across a wide range of scales. Hierarchi- 
cal clustering no longer provides a guiding principle for theoretical analyses and 
the computation time required to carry out credible simulations becomes pro- 
hibitively high. To gain insight into this problem, we perform high-resolution 
(N = 720 3 — 1584 3 ) simulations of an Einstein-de Sitter cosmology where the 
initial power spectrum is P(k) oc k n , with —2.5 < n < —1. Self-similar scaling is 
established for n = — 1 and n = —2 more convincingly than in previous, lower- 
resolution simulations and for the first time, self-similar scaling is established for 
an n = —2.25 simulation. However, finite box-size effects induce departures from 
self-similar scaling in our n = —2.5 simulation. We compare our results with 
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the predictions for the power spectrum from (one-loop) perturbation theory and 
demonstrate that the renormalization group approach suggested by McDonald 
(2007) improves perturbation theory's ability to predict the power spectrum in 
the quasilinear regime. In the nonlinear regime, our power spectra differ signif- 
icantly from the widely used fitting formulae of Peacock & Dodds (1996) and 
Smith et al. (2003) and a new fitting formula is presented. Implications of our 
results for the stable clustering hypothesis vs. halo model debate are discussed. 
Our power spectra are inconsistent with predictions of the stable clustering hy- 
pothesis in the high-A; limit and lend credence to the halo model. Nevertheless, 
the fitting formula advocated in this paper is purely empirical and not derived 
from a specific formulation of the halo model. 

Subject headings: Methods: N-body simulations — cosmology: dark matter 

1. Introduction 

In the standard cosmological paradigm, present-day structures arise from small-amplitude 
density perturbations which have their origin in the very early Universe. These primordial 
perturbations are presumed to form a Gaussian random field whose statistical properties are 
described entirely by the power spectrum, P{k). While higher order statistics are required 
to describe the density field once nonlinearities develop, the power spectrum remains central 
to our understanding of structure formation. 

During the radiation-dominated phase of the standard cold dark matter (CDM) scenario, 
the power spectrum evolves from its primordial, approximately power-law form, P(k) oc k 
to one in which its logarithmic slope, n e fj = d\nP(k)/d\iak, decreases from n eff ~ 1 on 
large scales to n c g ~ —3 on small scales 1 . At the start of the matter-dominated phase, 
which signals the beginning of structure formation, the dimensionless power spectrum, 
A 2 (A;) oc k 3 P(k), decreases monotonically with scale. The implication is that structure 
forms from the bottom up. Hierarchical clustering, as this process has come to be known, 
is the central idea in our understanding of structure formation. Hierarchical clustering also 
explains why cosmological N-body simulations are able to provide a reasonable facsimile of 
true cosmological evolution; with enough dynamic range, simulations are able to follow the 



1 More precisely, the logarithmic slope decreases from n c s — n p to n e ff = n p — 4 + logarithmic corrections 
where n p is the spectral index of the primordial power spectrum. An analysis of the three-year WMAP data 
by Spergel (2007) indicates that n p — 0.958 ± 0.016. However, for the sake of argument, we will set n p = 1 
since the difference is not relevant for our discussion. 
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development of virialized, highly nonlinear systems on small scales while properly model- 
ing the large-scale tidal fields that shape them. However, the dynamic-range requirement 
becomes increasingly difficult to achieve as n c g — > —3 since, in this limit, A 2 becomes in- 
dependent of k and structures collapse nearly simultaneously across a wide range of scales. 
Put another way, as n c fr — * —3, the infrared divergence of the power spectrum becomes 
increasingly problematic for numerical (as well as theoretical) studies. 

Interest in the small-scale limit of the CDM hierarchy was prompted by the realization 
that dark matter halos have a wealth of substructure (Moore, et al. 1999a; Klypin 1999) 
and that this substructure may have important implications for both direct and indirect 
dark matter detection experiments (See, for example, Stiff, Widrow, & Frieman (2001); 
Diemand, Kuhlen & Madau (2007); Kuhlen, Diemand, & Madau (2008) and Kamionkowski 
& Koushiappas (2008)). High-resolution simulations suggest that the subhalo mass function 
extends down to the dark matter free-streaming scale with approximately constant mass in 
substructure per logarithmic mass interval. These simulations probe structures which form 
from an initial power spectrum with —3 < n c fr < —2. For example, in the simulation of the 
first CDM objects by Diemand, Kuhlen & Madau (2007), n c g ~ —2.8. With such extreme 
spectra, care must be taken in order to insure that the results are not corrupted by finite- 
volume effects. This issue, as it relates to the halo and subhalo mass functions, is discussed 
in Power & Knebe (2006) and Bagla & Prasad (2006) as well as in the companion to this 
paper, Elahi et al. (2008). 

In this paper, we provide insight into the small-scale limit of CDM by focusing on scale- 
free cosmologies, that is Einstein-de Sitter cosmologies where the initial power spectrum is 
a power-law function of k, P{k) oc k n . The guiding principle for understanding structure 
formation in these models is self-similar scaling which implies that the functional form of 
the dimensionless power spectrum is time- independent, up to a rescaling of the wavenumber 
k (see below). Self-similar scaling provides a diagnostic test of whether a simulation has 
sufficient dynamic range (Jain & Bertschinger 1998). Contact with the standard ACDM 
cosmology is made by treating parameter n as a proxy for scale: n ~ —1.8 corresponds to 
cluster scales and n ~ —2.2 to galactic scales. The limit n — > — 3 corresponds to the bottom 
of the CDM hierarchy. 

We perform N-body simulations with n = — 1, —2, —2.25, and —2.5, and compare our 
results directly with theoretical models. The computation costs to conduct credible simula- 
tions with n < —2.5 are prohibitively high (see below) and even our n = —2.5 results must be 
considered suspect because of finite box-size effects. We compare our results for the power 
spectra in the quasilinear regime with predictions from one-loop perturbation theory and 
demonstrate that for n = — 2 and —2.25, the agreement is vastly improved if one implements 
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the renormalizat ion-group approach suggested by McDonald (2007). 

To predict the full non-linear power spectrum, one must resort to semi-analytic models 
such as the ones described in Hamilton et al. (1991) and Peacock & Dodds (1996). These 
models are based on the stable-clustering hypothesis (Peebles 1974; Davis & Peebles 1977) 
which holds that gravitationally-bound systems decouple from the rest of the Universe once 
they collapse. An alternative, known as the 'halo model' (Ma & Fry 2000; Peacock & Smith 
2000; Seljak 2000), allows for the continual accretion of mass onto existing haloes. The 
density field is treated as a distribution of mass concentrations, each characterized by a 
density profile. The power spectrum then involves the convolution of this density profile 
with the halo mass function. Simulations by Smith et al. (2003) of structure formation in 
a number of scale-free cosmologies demonstrate a clear departure from the stable clustering 
hypothesis and appear to support the halo model, at least qualitatively. 

Peacock & Dodds (1996) and Smith et al. (2003) provide fitting formulae for nonlinear 
power spectra. Formally, these formulae apply to all initial power spectra with n > — 3 
though they are calibrated using simulations with n > —2. One goal of this paper is to 
provide an alternative fitting formula which applies when n < — 2. 

The overall layout of this paper is as follows: In §2, we present background material 
including a discussion of self-similar scaling, perturbation theory, and semi-analytic models. 
We describe our simulations in §3 and our results in §4. We also provide a new and improved 
fitting formula for the nonlinear power spectra. In §5, we discuss the implications of our 
results for the halo model and stable clustering hypothesis. We conclude, in §6, with a 
summary and a discussion of directions for future investigations. 



In keeping with standard definitions (e.g. Peacock 1999), we express real-space density 
perturbations as deviations from the mean background density, Phg{t), and then construct 
the /c-space representation as follows: 



2. 



Preliminaries 



2.1. Statistics of the Density field 



<J(x, t) 
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The power spectrum, P(k), (strictly speaking a spectral density) is defined by the relation, 

(2nf5 D (k-k')P(k') = (5(k)5(k')), (3) 

where (• • • ) denotes an ensemble average, 5d is the Dirac delta function, and where the 
t-dependence is dropped for notational simplicity. Note that both 5 (k, t) and P (k, t) have 
units of volume. So long as the density field is statistically isotropic, k — \k\ encapsu- 
lates the wavenumber dependence. A useful quantity is the dimensionless power spectrum, 
A 2 (k) = k 3 P(k)/2ir 2 , which measures the power per logarithmic wavenumber bin. Given 
these definitions, the wavenumber at which A 2 ~ 1 corresponds to the dividing line between 
the linear and nonlinear scales. 



2.2. Initial Conditions 

Initial conditions for our N-body simulations are specified at an early epoch when the 
density field is accurately described by linear perturbation theory. We assume that prior to 
this epoch, the power spectrum is a power-law function of k between appropriately chosen 
high and low wavenumber cutoffs. That is, 

P (k, a) = P L (k, a) f uv (k, k c ) f IR (k, k B ) (4) 

where 

P L (k,a) = Aa 2 k n , (5) 

a is the scale factor, and A is a normalization constant. (Here and throughout, we assume 
an Einstein-de Sitter cosmology. For other cosmological models, a is replaced by the linear 
growth factor, D(a).) The functions fuv an d fm truncate the power spectrum above k c 
and below k B , respectively. We assume k B = 2n/B and k c = k^ y = nN 1 ^ /B where B is 
the size of the simulation "box" and k^ y is the Nyquist wavenumber of the initial particle 
distribution. We also follow Kudlicki, Plewa, & Rozyczka (1996) in using the truncation 
functions 

fuv(k, k c ) = e -( fe W f IR (k, k B ) = e(k - fc fl ). (6) 

where Q(x) is the Heaviside function. A similar initial set-up is used in our renormaliza- 
tion group calculations (see section 4.1). Our assumptions imply that at this early epoch 
A 2 (k c ) <^ 1, so that all modes within the simulation or computation box are initially in the 
linear regime. 

Utilizing these definitions, we describe the nonlinear scale L NL = 2n/kNL by the con- 
dition A 2 L (k = kNL, a) = 1, i-e-, k NL = (47rAa 2 )~ 1 ^ n+3 \ Since L NL is the only preferred 
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length scale, the power spectrum should evolve according to the self-similar scaling ansatz 

A 2 (k,a) = A 2 (k/k NL ) (7) 

(see Jain & Bertschinger (1996, 1998) for a discussion and review of the literature). Note 
that since A| = {k/k NL ) n+3 , we can also write A 2 as a function of A|. 

2.3. Perturbation Theory 

The evolution of the power spectrum from the linear regime to the mildly nonlin- 
ear or quasilinear regime can be estimated via perturbation theory (PT, see, for example, 
Bernardeau et al. (2002)). The starting point is an expansion for the density perturbation 
field of the form: 

oo 

5(k,t) = J2 = * n (t)5 n (k) (8) 

n=l 

where a is again the scale factor, 5\ characterizes linear density fluctuations, and 5 n denote 
terms of order (5i) n . So long as the density field is statistically isotropic, the perturbative 
expansion can be written in terms of the power spectrum: 

Ppt(M) =P L (k,a) + pV(k,a) + ... (9) 

where = O (A; 3 P|). In practice, the series is rarely carried beyond second order in P L . 
The calculations are aided by a diagrammatic scheme in which the higher-order terms are 
described as "loop corrections" to the "tree- level" term, P L (Scoccimarro & Frieman 1996). 
The one-loop correction, P^\ comprises two distinct terms or diagrams, P i3 and P 22 . These 
terms involve integrals over the linear power spectrum, Pj,. Explicit expressions can be found 
in Scoccimarro & Frieman (1996), McDonald (2007) and elsewhere. For scale-free models, 
they combine to give 

A 2 PT (k, a) = A 2 L (k, a) (l + X(n)A 2 L (k, a)) + O ((A 2 ,) 3 ) (10) 

where A(n), which can be calculated analytically, is positive for n > —1.4 and negative 
for n < —1.4 (Scoccimarro & Frieman 1996; Bernardeau et al. 2002). Hence, n ~ —1.4 
represents a "critical index" where nonlinear corrections are vanishingly small (Scoccimarro 
1997). 

2.4. Renormalization Group Approach 

PT breaks down when the loop corrections, which are formally divergent, become com- 
parable to the tree-level terms. The renormalization group (RG) scheme proposed by Mc- 
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Donald (2007) alleviates this problem essentially by updating the perturbative expansion 
as the system evolves. RG removes the divergences in the PT expansion, leaving behind a 
well-behaved, renormalized power spectrum. Operationally, one begins with an initial power 
spectrum and takes a small step forward in time using one-loop perturbation theory. The 
new power spectrum is used as the initial condition for the next timestep. This procedure 
is accomplished by solving the integro-differential equation 

dP r , ■ l'li (ii) 



d(a 2 ) 

with the initial conditions, P = Pl/o 2 - Here, P± 3 = Piz/a 2 and P 2 i = Pn/a 2 with the 
proviso that P rather than P L is used in evaluating P 13 and P 22 (see McDonald (2007)). 
Details on our own scheme for solving this equation are given in Section 4.1. 

As discussed in McDonald (2007), the method has a number of limitations. In particular, 
Eq. 11 ignores higher-order terms (two- loop and beyond) in P L . Moreover, decaying mode 
solutions are not included. Thus, our RG results should be interpreted with a degree of 
caution. 

The approach described here is an example of how RG techniques can be used to re- 
move secular divergences in differential equations. Crocce & Scoccimarro (2006) and Crocce 
& Scoccimarro (2008) outline an alternative method to study the nonlinear evolution of 
large-scale structure which also employs RG techniques. Their formalism is conveniently 
represented in terms of Feynmann diagrams and is closer in spirit to RG applications in high 
energy and statistical physics. 



2.5. Nonlinear Regime — Stable Clustering vs. Halo Model 

N-body simulations provide the most direct means to determine the nonlinear power 
spectrum. The so-called HKLM procedure provides an avenue by which one can begin to 
understand the simulation results from a theoretical basis (Hamilton et al. 1991; Peacock 
& Dodds 1996). The approach yields fitting formulae for the nonlinear power spectrum (or 
two-point correlation function, £(r)) in different cosmological models. The key assumption is 
the existence of a one-to-one relation between the nonlinear power spectrum at wavenumber 
k and the linear power spectrum at an earlier epoch and smaller wavenumber, k^. The 
starting point is the mapping 



k L =(l + A 2 (k,a)) k, 



(12) 
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which is the analogue of the real-space relation, 



r L = (l + Z(r)y'r, 



(13) 



that results from associating the volume averaged correlation function, £, with the local 
over- density. The nonlinear power spectrum may then be written as a function of the linear 
power spectrum at the earlier epoch: 



where / is an appropriately chosen fitting formula and, by assumption, A 2 (/c^) <C 1. Note 
that for scale-free models, Eq. 12 together with Eq. 7 implies Eq. 14. In the linear regime 
f(x) must be the identity function, i.e., f(x) = x, and the power spectrum grows with the 
scale factor as A| oc a 2 . In the nonlinear regime, Hamilton et al. (1991) and Peacock & 
Dodds (1996) appeal to the "stable clustering hypothesis", which holds that highly nonlinear 
structures decouple from the expansion. Under this assumption, A 2 oc a 3 and the asymptotic 
behaviour of f(x) must be given by / oc x 3 / 2 . 

To estimate f(x), Peacock & Dodds (1996) advocate the fitting formula 



where the parameters 7, (3, V, A, and B are determined by fitting Eq. 15 to power spec- 
tra measured in simulations. The parameters are understood as follows: B determines a 
second-order departure from linear growth, A and 7 control the behaviour in the quasilinear 
regime, V controls the amplitude of the asymptote and (5 shapes the transition between 
the two regimes. Best-fit parameters are expressed as functions of 1 + n/3, for example, 
7 = 3.310 (1 + n/3) -0 ' 244 . The expressions for the other parameters similarly diverge as 
n — > — 3 though one must bear in mind that they are based on results from simulations with 



The assumption of stable clustering has been challenged by various groups (Ma & Fry 
2000; Peacock & Smith 2000; Seljak 2000; Smith et al. 2003) on the basis that dark matter 
haloes continually accrete matter and never fully decouple from the rest of the Universe. 
An alternative approach is provided by the halo model in which the density field is given 
as a distribution of mass concentrations (haloes) which evolve and have their own internal 
structure. The two-point correlation function comprises a one-halo term, which is associated 
with the correlation of mass within a single halo, and a two-halo term, which is associated 
with the correlation between different haloes (Ma & Fry 2000; Peacock & Smith 2000; Seljak 
2000; Smith et al. 2003). Since the power spectrum is the Fourier transform of the two-point 



A 2 



(M) = / (A|(* L )) , 



(14) 




(15) 



n > -2. 
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correlation function, the associated components of the power spectrum, P\h and P2/1, involve 
integrals over the halo mass function and Fourier-transformed halo density profile. We return 
to this point in Section 5. 

Motivated in part by the separation of components used in the halo model, Smith et al. 
(2003) construct a fitting formula in the form 

A 2 (k, a) = Aq (k, a) + A 2 NL (k, a) (16) 

for the power spectra measured in their simulations. By construction, Aq dominates the 
power spectrum in the quasilinear regime and is meant to account for halo-halo correlations 
while A 2 NL dominates the power spectrum in the nonlinear regime and is meant to account for 
single halo correlations. However, the model is purely empirical and not calculated directly 
from the halo model. It is also worth noting that their formula has eight free parameters, 
three more than that of Peacock & Dodds (1996). 



3. Simulations 

N-body simulations are carried out with scale-free initial conditions and n = —1, —2, 
—2.25, —2.5 using the parallel tree-PM code GADGET-2 (Springel 2005). Initial conditions 
are generated on a regular grid using a second order-Lagrangian Perturbation Theory (2LPT) 
code (Crocce, Pueblas, & Scoccimarro 2006; Thacker & Couchman 2006). The primary 
benefit of 2LPT is to reduce the impact of spurious transient modes which arise from the 
truncation of the perturbative expansion. Since these modes decay, their impact is to delay 
the time in the simulation at which credible statistics can be calculated (Scoccimarro 1998; 
Crocce, Pueblas, & Scoccimarro 2006). The simulations are run with a softening length of 
1/30 the initial interparticle spacing. GADGET-2 parameters such as the opening angle 
used in constructing the particle tree and maximum time step criterion were set to their 
default value. 

Table 1 summarizes key features of the simulations used in this study such as the 
number of simulation particles and the epochs at which the power spectra are measured. 
The latter are expressed in terms of the ratio a/a* where a* is defined as the scale factor 
at the epoch when the mode on the scale of the box is equal to the nonlinear scale, that is, 
a/a, = (k B /k NL ) (n+3)/2 . With this definition, (a/a,) 2 = A 2 (k B ). 

As n approaches —3, the absence of modes beyond the box scale induces an error in the 
nonlinear power spectrum. Smith et al. (2003) adopt the ad hoc criterion that the missing 



variance, a m i SS , associated with these modes satisfies the condition 



miss — 



< 0.04 



(17) 



where, to a good approximation, a miss = A| (k B ) G(3 + n) with G(y) = (1 — 0.31y + 0.015|/ 2 + 
0.00133y 3 )/y. Figure 1 shows a miss for the outputs of our six simulations. We see that the 
final two outputs in our high-resolution n = —2 and n = —2.25 simulations and all but the 
first few outputs in our n = —2.5 simulation fail to meet the Smith et al. (2003) criterion. 
We return to this point below. 

Power spectra are calculated using a cubic mesh with side length L (Table 1). We set 
L = 2N rn so long as N m is a power of 2. Here, N m = N 1 ^ 3 is the side length of the initial 
grid of simulation particles. When N m is not a power of 2 we set L equal to the first power 
of 2 larger than N m . The power spectra are calculated using piecewise quadratic spline 
interpolation (Hockney & Eastwood 1981) and adjusted to account for the strong filtering 
of this mass-assignment scheme. No correction is made for shot noise. 

The ratio of the dimensionless power spectrum at the Nyquist frequency, k^ y , to that 
at the box scale, k B , provides a measure of a simulation's dynamic range. For a scale-free 
power spectrum 



Thus, if N = 256 3 particles are required to achieve a scaling solution over a reasonable range 
in A 2 when n = —2, (Jain & Bertschinger 1998), 256 4 ~ 1625 3 particles are required at 
n = —2.25, 256 6 particles are required for n = —2.5, and 256 12 particles are required for 
n = —2.75. We set the particle number for the n = —2.25 simulation on the basis of these 
arguments. 

We fully anticipated that self-similar scaling would not be achieved in our n = —2.5 
simulation. Our run, carried out with iV = 720 3 illustrates the difficulties that arise as one 
attempts to simulate highly negative spectral indices. In practice iV = 1584 3 is the largest 
simulation we can perform within word-addressing limits. It is also clear from this discussion 
that running an n = —2.5 simulation at 1584 3 will yield little improvement over our 720 3 
simulation since, apparently, one requires iV = 65536 s . Further discussions of the difficulties 
in simulating scale-free models with n — > —3 can be found in Smith et al. (2003) and Elahi 




(18) 



et al. (2008). 



4. Results 



In Figure 2 we plot the power spectrum from our N = 720 3 n — — 1 simulation at the six 
epochs listed in Table 1. The power at a given wavenumber increases with time while k^L 
(in this figure, the wavenumber where the power spectrum deviates from the linear form, 
P(k) oc A; -1 ), decreases. These results are in close agreement with previous simulations. The 
cumulative halo distribution is also in good agreement with the expected results (Elahi et 
al. 2008). In Figure 3 we test the self-similar scaling ansatz by plotting the dimensionless 
power spectra as a function of k/kNL for each of our four high-resolution simulations. Taken 
together, the spectra from different epochs yield a composite power spectrum. Consider, 
first, the case n = — 1 (upper left panel). The power spectrum covers 7 orders of magnitude 
in A 2 or, equivalently, 3-4 orders of magnitude in k. The fact that the composite power 
spectrum is very nearly a single- valued function of kjk^L indicates that self-similar scaling 
is essentially achieved. Note also that A 2 /A| < 1 for most values of k. This result is 
consistent with PT (See Eq. 10 and note that A(n = —1) < 0). 

As with the n — — 1 run, the composite spectra for n = —2 and n = —2.25 show excellent 
consistency with the scaling hypothesis. However, a departure from self-similar scaling is 
observed at large a in the n = —2.5 simulation highlighting the difficulty of simulating 
n — > — 3 spectral indices. Note that the high-A; feature in some of the early timesteps of 
our n = —2.5 simulation is a remnant of the grid used in setting up the initial conditions. 
The feature is subdominant to the physical small-scale power at later times in the n = —2.5 
simulation. 

The effect of resolution in achieving self-similar scaling is illustrated in Figure 4 where 
we compare the spectra from the three n = — 2 simulations. A departure from self-similar 
scaling is apparent in the N = 256 3 simulation and quite severe for N = 32 3 . It may be that 
these simulations are over-evolved, as suggested by Figure 1. In any case, based upon the 
scaling arguments presented in the previous section, an N = 32 3 , n = —2 simulation should 
be comparable to an N ~ 32 6 = 1024 3 , n = —2.5 simulation. Hence, it is not surprising 
that the departures from self-similar scaling seen in the right-hand panels of Figure 4 are 
comparable to those seen in our n = —2.5 simulation (lower-right panel of Figure 3. The 
departure manifests itself in a suppression of power at small k or large scales, as expected 
since power is missing due to the finite size of the simulation volume. 
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4.1. Perturbative or Quasilinear Regime 

We now focus on the quasilinear regime in order to illustrate the improvement renor- 
malization group methods brings to perturbation theory. We implement the RG approach 
by solving the Eq. 11 assuming an initial scale-free power spectrum with n = —1, —2, —2.25 
or —2.5. The initial spectrum is "evolved" forward in time (or equivalently, scale factor a) 
using a 4th-order Runge Kutta scheme with an adaptive stepsize (see, for example, Press 
et al. (1992)). Each Runge Kutta step requires an evaluation of the Pi 3 and P 2 2 integrals. 
These integrals have the same functional form as those that appear in Scoccimarro & Frie- 
man (1996) except that Pl is replaced by P which, in turn, is updated at each step. As 
with N-body simulations, we must truncate the initial power spectrum at both high and 
low wavenumbers. Otherwise, the integrals would diverge. Scoccimarro & Frieman (1996) 
use sharp cutoffs which are convenient for power-law spectra with integer n since analytic 
expressions can be derived. Smooth cutoffs are more manageable for the RG analysis where 
P 13 and P 22 must be evaluated numerically (Scoccimarro, private communication) . 

To make contact with our discussion in Section 2.2 we refer to the IR and UV cutoffs 
as ks and k c , respectively and assume an initial power spectrum of the form 

P L (k) = Aa 2 k n f uv (k, k c , A) f IR (k, k B , A) (19) 

where 

/„ (*, k B , Ak) = i (erf (!^>) + l) fuv (*, K A*) = \ eric (!^) . 

(20) 

The ratio k c /ks, which corresponds to the dynamic range of the calculation, is set to 10 6 
while A, which determines the sharpness of the k-space cutoffs, is set equal to 0.5. For 
n < —1, P13 and P 2 2 have terms of order (/ce/A; c ) n+1 which cancel, leaving behind a residual 
term of order k 2n+3 . The challenge, numerically, is to determine the surviving terms which 
can be much smaller than the terms that cancel, especially for large k and small n. Here, 
we use the Romberg integration routine from Press et al. (1992). 

The solution to Eq. 11 yields an evolutionary sequence for the power spectrum, P(k, a). 
Departures from the linear power spectrum increase with a beginning at high wavenumber. 
We evolve P until the k^L is roughly equal to the geometric mean of k c and ks- Since 
our dynamic range is a full three orders of magnitude greater than is found in our N- 
body simulations, finite box effects are much less a concern here. Moreover, our results are 
insensitive to the form of the cutoff functions, fuv an d fm-, since they are derived in a region 
well inside the computation box. 

In Figure 5 we show the measured A 2 in the mildly nonlinear regime together with 
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predictions from PT, RG- improved PT, and the Zel'dovich approximation. The latter is 
discussed in Taylor & Hamilton (1996). Also shown are the fitting formulae of Peacock 
& Dodds (1996) and Smith et al. (2003). Note that in the n = -1 case, A 2 /A 2 NL slowly 
decreases with increasing k for k < k^L and it is difficult to discern the true self-similar 
evolution of the power spectrum given that the actual initial power spectrum has a large-A; 
cutoff. On the other hand, for n = —2, and —2.25, it is clear that RG does the best job 
of tracking the power spectrum in the quasilinear regime the RG power spectra does not 
quite capture the rapid evolution of the measured power spectra. The situation is less clear 
for n = —2.5 where the validity of the simulation is in doubt. The question remains as to 
whether agreement between RG and the simulations might be improved by refinements in 
the RG analysis, such as those suggested by McDonald (2007). 



4.2. Nonlinear Regime 

In Figure 6, we plot the complete power spectrum for our high-resolution, n — — 1 
simulation together with the predictions of Peacock & Dodds (1996) and Smith et al. (2003). 
Also shown is our own fitting formula given by 

A 2 (A;) = A|(%(fc/fc NL ) (21) 

where 

f l + Ax + Bx a Y 

9{X)= { 1 + ) ■ (22) 

The form of this formula is motivated by that of Peacock & Dodds (1996) but has one 
additional parameter to allow for a more general behaviour in the k — > oo limit. The 
parameters, derived by performing a nonlinear least-squares fit (see, for example, Press et 
al. (1992)), are given in Table 2. The lower panel in Figure 6 shows the logarithmic slope 
pi = d\nP/d\nk = d\nA 2 /d\nk — 3. Note that \x monotonically decreases with k near the 
Nyquist wavenumber. 

In Figures 7 and 8, we show A 2 and /i for our high-resolution, n = — 2 and —2.25 
simulations. Again, \i decreases monotonically in the large- A; limit. Nevertheless, by design, 
virtually all fitting formulae (including our own) have a power-law form in the high- A; limit, 
that is, P(k) oc k^ as k — > oo. The lesson is that fitting formula should not be extrapolated 
to scales below the smallest scales probed by the simulation used in their construction. 

Neither the Peacock & Dodds (1996) nor Smith et al. (2003) fitting formulae do a 
particularly good job of fitting the power spectra from our simulations. For n = —2, the 
Smith et al. (2003) formula provides a reasonable fit up to k/k NL ~ 5 but decreases too 
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rapidly beyond this point. Conversely, Peacock & Dodds (1996) predict that A 2 /A| is 
constant in the large-A; limit whereas the measured power spectrum shows a clear decline. 
The discrepancies between predicted and measured power spectra for n = —2.25 are equally 
severe. By contrast, Eq. 21 provides an excellent fit to the nonlinear power spectra from our 
high-resolution simulations. And while it has one more parameter than the fitting formula 
of Peacock & Dodds (1996), it has two fewer than that of Smith et al. (2003). 

The n = —2.5 case, shown in Figure 9, is difficult to analyse because of the departure 
from self-similar scaling. In this plot, we truncate the power spectra from different timesteps 
at large k, where the effects of aliasing is apparent, and at small k, where the effects of 
missing large-scale power is apparent. We contend that this procedure yields a roughly 
continuous curve, which provides a reasonable facsimile of the true power spectrum. The 
plausibility of this procedure is illustrated in the left-hand panels of Figure 4 where one can 
imagine carrying out a similar procedure with our n = — 2, N = 256 3 results to yield an 
approximate form for the power spectrum from our high-resolution simulation. 



5. Halo Model Revisited 

In this section, we explore the halo model vis-a-vis our simulation results in more detail. 
Our focus here is on the high-A; limit of the nonlinear power spectrum where the one-halo 
term dominates. The term can be expressed as an integral over the halo mass function, 
dn/dM, and the Fourier-transformed density profile of a single halo, p(k, M). Following 
Seljak (2000), we use the peak height v = (5 C / o (M)) 2 as the integration variable where 5 C 
is the critical overdensity for spherical collapse (S c ~ 1.68 in an Einstein-de Sitter universe) 
and <t(M) is the rms mass overdensity for a spherical region of radius R = (3M/47rpb g ) 1 ^ 3 - 
For scale- free models, M oc z/ 3 /( n+3 ). One finds 

, 1 2 



(2*) 

where 



(-) 


-~p(k, M)l 


VPbg/ 


M 



dv . (23) 



, . , M dn dM /ni . 

ftM = ^w (24) 

is a dimensionless form for the halo mass function. 

A commonly used fitting formulae for halo density profiles take the form 

P( r ) = , , x 7n P " , ( 25 ) 

{r/r s ) (1 + r/r s ) 

where 77 ~ 3 and 7 ~ 0.5 — 1.5 (Navarro, Frenk, & White 1996; Moore, et al. 1999b; Kravtsov 
ct al. 1998). The characteristic halo scale length, r s , depends on the halo mass through 
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Table 1. Summary of Simulations 



n N L initial scale factor (aj /a*) output scale factors (a/a*) 



-1 720 3 1024 0.0014 

-2 32 3 64 0.028 

256 3 512 0.010 

1024 3 2048 0.006 

-2.25 1584 3 2048 0.009 

-2.5 720 3 1024 0.015 



0.026, 0.11, 0.19, 0.21, 0.24, 0.30 
0.052, 0.13, 0.23, 0.31, 0.42, 0.49 
0.031, 0.054, 0.17, 0.29, 0.51, 0.67 
0.024, 0.043, 0.063, 0.17, 0.30, 0.36 
0.018, 0.021, 0.039, 0.084, 0.221,0.394 
0.033, 0.074, 0.15, 0.23, 0.26, 0.29, 0.33 



Table 2. Parameters for fitting formula, Eq. 21 



A B C a 7 (3 



-1 -0.158 0.181 0.0729 1.571 1.845 2.664 

-2 -0.0312 0.690 0.478 1.243 1.266 8.647 

-2.25 3.471 4.038 0.348 1.413 1.372 0.659 

-2.5 174.0 110.2 4.532 1.492 1.231 0.879 
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the relation r s = r vir /c vir where r vir is the virial radius and c v i r is the halo concentration 
parameter. Cosmological simulations indicate that the concentration parameter decreases 
with mass, roughly as a power-law (see, for example, Navarro, Frenk, & White (1996) and 
Bullock et al. (2001)). With this in mind, we follow Seljak (2000) and Ma & Fry (2000) and 
write c vir = (M vir /Mo) _/3 so that r s oc M^ 3 ^' 3 . 

Our focus is on the power spectrum in the high-A; limit where the Fourier transform of 
p(r) may be approximated by a step function, 

p(k, M) ~ MO (1 — kr s (M)) . (26) 

(At wavenumbers above k = r^ 1 , p decreases as (Ax s ) 3-7 but the Heaviside function provides 
a suitable form for our discussion.) Press & Schechter (1974) and Sheth & Tormen (1999) 
provide analytic expressions for h(u). In the small- v limit (i.e., small M or large k limit), 
one finds vh oc v a where a = 0.5 (a = 0.3) for Press & Schechter (1974) (Sheth & Tormen 
(1999)). Putting all this together, we find that in the large-/c limit Aj h oc k^ lh+3 where 

Hu> ~ ^ 3 (27) 

Ma & Fry (2000). Note, however, that pih is independent of the cusp parameter 7 indicating 
that the power spectrum in the strongly nonlinear regime is insensitive to the structure of 
the inner halo. 

In Figure 10 we plot our results for the asymptotic slope of the nonlinear power spectrum, 
p, together with those from the Smith et al. (2003) simulations. We make the point of 
displaying our results as upper bounds on p since the slope of p appears to be a decreasing 
function of k as k — > k^ y (see Figures 6-9). For n — — 1 and —2, these bounds agree with 
the quoted values from the Smith et al. (2003) simulations. Furthermore, the functional 
dependence of p on n from their fitting formula appears to be consistent with our n = —2.25 
and n = —2.5 results. 

Our results suggest that p increases with increasing n and that p(n — > —3) = —3. In 
other words, as n — > —3, the nonlinear dimensionless power spectrum becomes indepen- 
dent of k (i.e., equal power per logarithmic wavenumber bin) just as with the linear power 
spectrum. Formulations of the halo model with constant, nonzero (3 cannot reproduce this 
behaviour. To illustrate this point, we plot these predictions for pih assuming (5 = 0.15, as 
in Seljak (2000), and either the Press & Schechter (1974) or Sheth & Tormen (1999) values 
for a. These predictions are inconsistent with our simulation results and those of Smith et 
al. (2003). 

Clearly, the dependence of the concentration parameter on halo mass is central to the 
development of the halo model. Bullock et al. (2001) devised a toy model to explain the 
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concentration-mass relation seen in simulations. In the case of scale-free cosmologies their 
model predicts (5 = (n + 3)/6 which, when combined with Eq. 27, would seem to yield the 
desired behaviour for p^h in the n — > —3 limit. The lower panel in Figure 10 shows this 
prediction. While it does better than constant-/? versions of the halo model, the predicted 
\x tends to lie above the values obtained in the simulation. 

6. Summary and Conclusions 

The fundamental tenet of the hierarchical clustering scenario is that small-scale objects 
form earlier than large-scale ones. A corollary of this statement is that individual structures 
can be identified with a specific wavenumber range of the primordial power spectrum ac- 
cording to their mass. In CDM cosmologies, the logarithmic slope or spectral index of the 
primordial power spectrum runs from 1 at large scales to -3 at small scales. Thus, as we 
increase the dynamic range in our simulations and push to smaller and smaller scales, we 
probe structures that form from density perturbations with a power spectrum approaching 
k~ 3 . However this limit represents a singular case where the dimensionless power spectrum 
is independent of scale and structures across a wide range in mass collapse nearly simul- 
taneously. The nature of structure formation changes and the computing requirements for 
performing simulations increase dramatically. 

This work and our companion paper, Elahi et al. (2008), provide inside into the under- 
lying physics of ACDM models by considering scale-free cosmologies. We focus here on the 
nonlinear power spectrum and in Elahi et al. (2008), on the distribution of subhaloes. The 
evolution of the power spectrum in scale-free cosmologies is remarkably simple — the dimen- 
sionless power spectrum, when written as a function of the ratio k/k^i, is time-independent. 
Obviously in simulations, the finite computation volume breaks the scale-free nature of the 
problem and leads to departures from the scaling solution. The dimensionless power spec- 
trum provides a simple test of whether finite volume effects have corrupted the simulation 
(Jain & Bertschinger 1998). 

Our high-resolution n = — 1, —2 and —2.25 simulations demonstrate the scaling solution 
across the simulation volume while showing clear differences with simulations performed at 
lower resolution. Moreover, our results differ markedly from the the fitting formula provided 
by Peacock & Dodds (1996) and Smith et al. (2003). A plausible power spectrum for n = 
—2.5 was constructed by stitching together outputs from different timesteps. Though it 
shows significant, and entirely expected departures from the scaling solution, it represents 
our best estimate of the power spectrum for models with n this small. We summarize our 
results for our four high-resolution simulations by means of a simple fitting formula. Future 
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work will fill in the gaps in n (e.g., n = —1.25, —1.5, —1.75, and —2.75) and enable use to 
develop a model for power spectra of arbitrary n and therefore arbitrary shape. 

The renormalization group improvements to perturbation theory developed in McDon- 
ald (2007) represent a promising avenue for studying scale-free models with n < — 2 and 
likewise, the low-mass limit of the CDM hierarchy. Not surprisingly, the calculations be- 
come more difficult as n — > —3. Our analysis of the n = —2 and —2.25 cases confirms 
McDonald's claim that RG does lead to an improvement in the predictions of perturbation 
theory when compared to simulations with the caveat that, as n becomes more negative, the 
RG-predicted power spectrum fails to capture the rapid rise of the power spectrum seen in 
the simulations. Our analysis for the n = —2.5 case is less conclusive but departures in the 
simulations from the scaling solution suggest that the problem may reside there rather than 
in the PT analysis. In principle, RG-improved PT can yield a handle on the form of the 
power spectrum in the mildly nonlinear regime even asn-> —3. 

The halo model provides a theoretical framework for understanding the two-point cor- 
relation function and nonlinear power spectrum. Our results, with respect to this model, 
are somewhat inconclusive. We agree with Smith et al. (2003) that the stable clustering 
hypothesis of Hamilton et al. (1991) and Peacock & Dodds (1996) fails. On the other hand, 
the halo model appears to have difficulty reproducing the relation between the asymptotic 
slope of the power spectrum and n. We must therefore settle for an empirical fitting formula 
for the nonlinear power spectrum. 
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Fig. 1. — Missing variance, cr miss , as given in Eq.17, for the six simulations described in 
the text with line types/colours as follows: solid/blue — n = — 1; dotted/red — n — —2, 
N = 32 3 ; short-dashed/red — n = -2, N = 256 3 ; long-dashed/red — n— —2, N = 1024 3 ; 
short-dashed-dot/green — n = —2.25; long-dashed-dot/cyan — n = —2.5. The output 
number, A^ output corresponds to the values listed in Table 1. For n = —2.5, we show cr miss 
for the last six outputs. The horizontal black curve corresponds to the criterion adopted by 
Smith et al. (2003). 
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Fig. 2. — Power spectrum, P(k) as a function of wavenumber k for the n = — 1 simulation. 
Different colours and symbols correspond to different output times as given in Table 1. The 
sequence, from early to late times is magenta-blue-cyan-green-brown-red, or, alternatively, 
open square-filled square-open triangle-filled triangle-open circle-filled circle. 
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Fig. 3. — Dimensionless power spectrum, A 2 as a function of wavenumber k. Colours and 
symbols are the same as in Figure 2. Bottom panel in each quadrant shows the ratio A 2 /A 2 . 
Upper left — n = — 1; Upper right — n — —2; Lower left — n = —2.25; Lower right - 
n = -2.5. 
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Fig. 4. — Dimensionless power spectrum, A 2 , and the ratio A 2 /A 2 for n = —2 from sim- 
ulations with different numbers of particles. Black points are power spectra at different 
timesteps measured in our highest resolution (N = 1024 3 ) simulation. Superimposed in 
colour are measurements from the iV = 256 3 (left) and N = 32 3 (right) simulations. As in 
Figure 2, the sequence, from early to late times is magenta-blue-cyan-green-brown-red or, 
alternatively, open square-filled square, open triangle-filled triangle-open circle-filled circle. 
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Fig. 5.— The ratio A 2 /A L clS £L function of k/k^i from our four high-resolution simulations 
as labelled in each panel. Black points are from the simulation with the different symbols 
representing results from different outputs (from early to late outputs: solid squares, open 
triangles, solid triangles, open circles, solid circles). Line colours/types are: blue/solid - 
RG; green/dotted — one-loop; cyan/dot-dashed — Zel'dovich approximation; red/short- 
dashed — Smith et al. (2003) fitting formula; magenta/long-dashed — Peacock & Dodds 
(1996) fitting formula. 
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Fig. 6. — The ratio A 2 /A^ as a function of k/k^L for the n = — 1, N m = 1024 simulation. 
Shown is the full range in k probed by the simulation is shown. Black points are from 
the simulation. Red curve is the Smith et al. (2003) fitting formula. Magenta curve is the 
Peacock and Dodds fitting formula. Blue curve is our own fitting formula, Eq. 21. Plotted 
in the lower panel is the logarithmic slope, /i of the power spectrum (see text). 
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Fig. 8. — Same as Figure 6 but for n = 



-2.25. 
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Fig. 9. — Same as Figure 6 but for n = 



-2.5. 
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Fig. 10. — Asymptotic slope of the power spectrum, ft, as a function of the slope, n, of the 
linear power spectrum. The black points are from our simulations while the green points are 
from the Smith et al. (2003) simulations. In the upper left panel we show the predictions of 
the halo model (Eq. 27) for f3 = 0.15 and a = 0.5 (Press & Schechter). The upper right panel 
shows the prediction of the halo model assuming /? = 0.15 and a = 0.2 (Sheth & Tormen). 
The lower left panel shows the predictions from Peacock & Dodds (1996) (magenta curve) 
and Smith et al. (2003) (red curve). The lower right panel, shows the predictions of the halo 
model using the prescription from Bullock et al. (2001) for (3. 



